Having visualized geodesics on surfaces and in space in previous presentations, it is high time to explore the gravitational field of a black hole. While this has been done before, generally for restricted cases on dead paper, the topic will benefit from the possibilities of live interactive visualization.

Begin with the Schwarzschild metric, which with the simplification $c=G=1$ has the form

$d{s}^2=(1-\frac{2M}{r})d{t}^2-(1-\frac{2M}{r}{)}^{-1}d{r}^2-r^{2}d{\theta }^2-r^2{sin}^2\theta d{\phi }^2$

This can be incorporated into a Lagrangian,

$L=\frac{1}{2}(\frac{ds}{d\tau }{)}^2=\frac{1}{2}[(1-\frac{2M}{r}){\stackrel{·}{t}}^2-(1-\frac{2M}{r}{)}^{-1}{\stackrel{·}{r}}^2-r^2{\stackrel{·}{\theta }}^2-r^2{sin}^2\theta {\stackrel{·}{\phi }}^2]$

from which one can evaluate geodesic equations using the Lagrange equations

$\frac{d}{ds}\frac{\partial L}{\partial {\stackrel{·}{q}}_i}=\frac{\partial L}{\partial q_i}$

The explicit equations are

$\begin{array}{c}\frac{d}{ds}\left[\right(1-\frac{2M}{r}\left)\stackrel{·}{t}\right]=0\\ \\ \frac{d}{ds}\left[\right(1-\frac{2M}{r}{)}^{-1}\stackrel{·}{r}]=-\frac{M}{r^2}{\stackrel{·}{t}}^2-\frac{M}{r^2}(1-\frac{2M}{r}{)}^{-2}{\stackrel{·}{r}}^2+r{\stackrel{·}{\theta }}^2+r{sin}^2\theta {\stackrel{·}{\phi }}^2\\ \\ \frac{d}{ds}\left(r^2\stackrel{·}{\theta }\right)=r^{2}sin\theta cos\theta {\stackrel{·}{\phi }}^2\frac{d}{ds}\left(r^2{sin}^2\theta \stackrel{·}{\phi }\right)=0\end{array}$

As with the previous presentations, one can expand these equations into second-order forms suitable for numerical integration,

$\begin{array}{l}\stackrel{··}{t}=-\frac{2M}{r^2}(1-\frac{2M}{r}{)}^{-1}\stackrel{·}{t}\stackrel{·}{r}\\ \stackrel{··}{r}=-\frac{M}{r^2}(1-\frac{2M}{r}){\stackrel{·}{t}}^2+\frac{M}{r^2}(1-\frac{2M}{r}{)}^{-1}{\stackrel{·}{r}}^2+(1-\frac{2M}{r})(r{\stackrel{·}{\theta }}^2+r{sin}^2\theta {\stackrel{·}{\phi }}^2)\\ \stackrel{··}{\theta }=-\frac{2}{r}\stackrel{·}{r}\stackrel{·}{\theta }+sin\theta cos\theta {\stackrel{·}{\phi }}^2\\ \stackrel{··}{\phi }=-\frac{2}{r}\stackrel{·}{r}\stackrel{·}{\phi }-2cot\theta \stackrel{·}{\theta }\stackrel{·}{\phi }\end{array}$

but since some of these equations are trivial to integrate, while at the same time providing physical context for initial conditions, it is appropriate to consider solutions that reduce the order of numerical integration.

The first equation allows one to introduce a constant of the motion, which being cognate to the temporal variable is an energy for the geodesic:

$(1-\frac{2M}{r})\stackrel{·}{t}=E\to \stackrel{·}{t}=E(1-\frac{2M}{r}{)}^{-1}$

Likewise the fourth equation alows the introduction of a second constant corresponding to a partial angular momentum for the geodesic:

$r^2{sin}^2\theta \stackrel{·}{\phi }=L_z\to \stackrel{·}{\phi }=\frac{L_z}{r^2{sin}^2\theta }$

Inserting this result into the third equation and then applying an integrating factor of $r^2\stackrel{·}{\theta }$ , one has

$\begin{array}{c}\frac{d}{ds}\left(r^2\stackrel{·}{\theta }\right)=\frac{L_z^2}{r^2}\frac{cos\theta }{{sin}^3\theta }\\ \\ \frac{d}{ds}[(r^2\stackrel{·}{\theta }{)}^2+\frac{L_z^2}{{sin}^2\theta }]=0\\ \\ (r^2\stackrel{·}{\theta }{)}^2+\frac{L_z^2}{{sin}^2\theta }=L^2\to \stackrel{·}{\theta }=\frac{1}{r^2}\sqrt{L^2-\frac{L_z^2}{{sin}^2\theta }}\end{array}$

where now a constant corresponding to total angular momentum for the geodesic appears.

These last two constants, as well as the differential equations that give rise to them, are identical in form to those appearing in the context of angular momentum in spherical coordinates. That means one can immediately write a conserved vector for motion in the presence of the black hole

$\mathbf{L}=r^2(\stackrel{·}{\theta }\hat{\mathbf{\phi }}-\stackrel{·}{\phi }sin\theta \hat{\mathbf{\theta }})$

that defines an invariant plane for geodesics with the given constants. This behavior is precisely that of the classical two-body Kepler problem due to the spherical symmetry of the system. This is also why most presentations of these geodesics use the restriction $\theta =\frac{\pi }{2}$ to simplify equations, but this is clearly not necessary. It is also inconvenient when visualizing geodesics in three dimensions: one wants to see the full behavior for various initial conditions.

The three first-order equations above can now be used to simplify the remaining Lagrange equation:

$\begin{array}{c}\frac{d}{ds}\left[\right(1-\frac{2M}{r}{)}^{-1}\stackrel{·}{r}]=-\frac{M}{r^2}(1-\frac{2M}{r}{)}^{-2}{E}^2-\frac{M}{r^2}(1-\frac{2M}{r}{)}^{-2}{\stackrel{·}{r}}^2+\frac{L^2}{r^3}\\ (1-\frac{2M}{r}{)}^{-1/2}\frac{d}{ds}\left[\right(1-\frac{2M}{r}{)}^{-1/2}\stackrel{·}{r}]=-\frac{M}{r^2}(1-\frac{2M}{r}{)}^{-2}{E}^2+\frac{L^2}{r^3}\\ \frac{d}{ds}\left[\right(1-\frac{2M}{r}{)}^{-1}{\stackrel{·}{r}}^2]=\frac{d}{ds}[(1-\frac{2M}{r}{)}^{-1}{E}^2-\frac{L^2}{r^2}]\end{array}$

using an integrating factor of $2\stackrel{·}{r}$ in the last step. There is now an additional constant of the motion,

$(1-\frac{2M}{r}{)}^{-1}{E}^2-(1-\frac{2M}{r}{)}^{-1}{\stackrel{·}{r}}^2-\frac{L^2}{r^2}=C$

which is proportional to the value of the Lagrangian and determines the nature of the geodesics. The constancy of the Lagrangian is to be expected, since having no ‘potential energy’ term it will be equal to the corresponding Hamiltonian, which is by definition constant.

Since this new constant is proportional to the original metric, the geodesics are time-like, null or space-like depending on whether $C>0$ , $C=0$ or $C<0$ . Adjusting the parametrization of the interval, this constant can be reduced to the values $C=\pm 1$ and $C=0$ .

In retrospect it is clearly true that the interval itself always provides a constant of the motion. For some reason this simple fact is generally not emphasized in textbooks...

The radial geodesic equation can now be written

$\stackrel{·}{r}=\sqrt{E^2-(1-\frac{2M}{r})(\frac{L^2}{r^2}+C)}$

The entire set of equations has been reduced in principle to a first-order system. In practice this is not entirely useful, since the presence of square roots will introduce branch cuts and hence spurious singularities into the solutions, creating difficulties in numerical integration. This is clearly to be avoided.

Combining all of these separate results, the system of equations to be integrated is

$\begin{array}{l}\stackrel{·}{t}=E(1-\frac{2M}{r}{)}^{-1}\\ \stackrel{··}{r}=-\frac{MC}{r^2}+\frac{L^2}{r^3}-\frac{3M{L}^2}{r^4}\\ \stackrel{··}{\theta }=-\frac{2}{r}\stackrel{·}{r}\stackrel{·}{\theta }+\frac{L_z^2}{r^4}\frac{cos\theta }{{sin}^3\theta }\\ \stackrel{·}{\phi }=\frac{L_z}{r^2{sin}^2\theta }\end{array}$

Note that with the introduction of the energy constant the explicit temporal variable falls out of the right-hand sides of equations. This variable, as a function of the parameter of interval along the geodesic, basically just reparametrizes motion along the geodesic. This behavior will only appear when visualizing geodesic motion temporally.

The parameter of interval and the temporal variable will both be naturally assumed to start at zero when evaluating geodesics numerically.

For other dynamic variables, it is desirable to specify a starting location in space of particular interest, just as one would for classical planetary orbits. Since the motion occurs in an invariant plane, one is tempted to infer the inclination of this plane to the xy-plane from ${\theta }_0$ , set ${\stackrel{·}{\theta }}_0=0$ at a turning point, and determine the partial angular momentum $L_z$ accordingly. This is a mistake, since it misses a large number of possible orbit orientations beyond those selected this way.

An arbitrary initial angle θ for the geodesic does not uniquely determine the inclination of the angular momentum vector to the z-axis, since the vector is free to rotate in space about the initial radial vector. The partial angular momement $L_z$ will have a simple relation to the angle θ at some point on the orbit, but not neccesarily the chosen initial point.

The correct procedure is to remember the full expression for the variation of this angular variable, namely

$\stackrel{·}{\theta }=\pm \frac{1}{r^2}\sqrt{L^2-\frac{L_z^2}{{sin}^2\theta }}$

The partial angular momentum must be left here as a free parameter, subject to the constraint of keeping this expression real.

The value of ${\phi }_0$ then allows rotation of the tilted invariant plane around the z-axis.

That just leaves initial conditions for the radial variable to be determined. As a consequence of the introduction of the remaining constants of the motion, the radial equation depends only on one variable. It is also quite interesting that this equation is no longer nonlinear for the same reason, now bearing a distinct similarity to classical motion. This does not get mentioned in textbooks...

Indeed, starting from a classical Kepler Hamiltonian with unit mass and deliberately mixed dynamic variables,

$H=\frac{1}{2}{\stackrel{·}{r}}^2+\frac{L^2}{2r^2}-\frac{M}{r}$

the corresponding radial equation of motion is

$\stackrel{··}{r}=-\frac{M}{r^2}+\frac{L^2}{r^3}$

For time-like geodesics, which are the ones physically realizable for massive bodies, the right-hand side of the general relativistic radial equation begins with precisely these two terms, then adds an inverse quartic correction. That means the general relativistic radial potential for time-like geodesics is

$V_{\mathrm{GR}}\left(r\right)=-\frac{M}{r}+\frac{L^2}{2r^2}-\frac{M{L}^2}{r^3}$

In order to have geodesics that circle the black hole, there must be a potential well just as for the classical case. Taking a derivative of the general relativistic potential, its extrema are located at

${V_{\mathrm{GR}}}^{\prime }=0\to r_{\mathrm{ext}}=\frac{L^2}{2M}(1\pm \sqrt{1-\frac{12M^2}{L^2}})$

In order for these extrema to remain real, the squared angular momentum must have the critical lower limit

$L_c^2\ge 12M^2$

Geodesics with squared angular momentum less than this value will decay into the singularity.

The value of the radial variable at the lower limit of the squared angular momentum is $r_c=6M$ . As the squared angular momentum increases one can expand the squared root term for the approximations

$r_{\mathrm{ext+}}\approx \frac{L^2}{2M}-3M-\frac{9M^3}{L^2}{r}_{\mathrm{ext-}}\approx 3M+\frac{9M^3}{L^2}$

The first term in the left-hand approximation is the result for classical Kepler orbits, providing circular orbits of arbitrary large radius both classically and general relativistically.

As the squared angular momentum grows without limit, only one extremum remains finite at $r_{\infty }=3M$ . Mathematically this is the same value as the single extremum for null geodesics, but curious all the same.

It will be useful to know the value of the general relativistic potential at the radial extrema. Rationalizing the value

$\frac{1}{r_{\mathrm{ext}}}=\frac{1}{6M}(1\mp \sqrt{1-\frac{12M^2}{L^2}})$

and performing a bit of algebra gives

$V_{\mathrm{GR}}\left(r_{\mathrm{ext}}\right)=-\frac{1}{6}+\frac{L^2}{108M^2}\pm (\frac{1}{9}-\frac{L^2}{108M^2})\sqrt{1-\frac{12M^2}{L^2}}$

A simple trick to help in this evalution is to use the intermediate substitution $x=M^2/L^2$ while remembering that it is a combination of squares of dynamic variables...

For motion that starts at a stationary radial position, the linear radial equation provides a relation between squared energy and other quantities:

${\stackrel{·}{r}}_0=0\to E^2=(1-\frac{2M}{r_0})(\frac{L^2}{r_0^2}+1)=1+2V_{\mathrm{GR}}\left(r_0\right)$

This general relativistic energy condition differs from that expected for a classical potential energy, but holds in this form nonetheless, and uniquely determines the arbitrary additive constant inherent in all potentials.

At the lower limit of squared angular momentum, the general relativistic potential and squared energy have the critical values

$V_{\mathrm{GR}}\left(r_c\right)=-\frac{1}{18}{E}_c^2=1+2V_{\mathrm{GR}}\left(r_c\right)=\frac{8}{9}$

which provides an absolute lower limit on energy for bound geodesics. An absolute upper limit is found by noting that $V_{\mathrm{GR}}\left(r\right)\to 0$ as $r\to \infty$ . In order for two turning points to be possible, the absolute maximum value of energy for bound geodesics is clearly unity.

For particular geodesics, the possible values of squared energy vary between the values at the two extrema, which are

$E^2\left(r_{\mathrm{ext}}\right)=\frac{2}{3}+\frac{L^2}{54M^2}\pm (\frac{2}{9}-\frac{L^2}{54M^2})\sqrt{1-\frac{12M^2}{L^2}}$

Following the expression for the lower limit of squared angular momentum, if this quantity is always denumerated proportional to squared gravitational mass, $L^2=\alpha M^2$ , then the practical limits for squared energy are

$E^2\left(r_{\mathrm{ext}}\right)=\frac{1}{54}(36+\alpha \mp \sqrt{\frac{(\alpha -12{)}^3}{\alpha }})$

For the lower limit $L^2=12M^2$ of squared angular momentum, this result is consistent with that above. With a bit of algebra, or just playing with numbers, the value for the lower sign here, corresponding to the smaller radial extremum, is unity when $L^2=16M^2$ . For larger angular momentum the value of squared energy remains capped by its absolute maximum of unity as stated above.

The constant M of course remains a free variables, setting the arena in which motion occurs, which is independent of particular geodesic motion. In practice the interactive graphics to follow are highly sensitive to the value of this variable, so it will be taken equal to one, $M=1$ , for simplicity of numerical evaluation and presentation.

All of these radial relationships can be summarized in a simple visualization of the general relativistic energy condition:

When visualizing radial motion it will make the most sense to start from the outer turning point and integrate inwards. This will ensure a connection to classical orbits, which exist for any initial radial value. The general relativistic geodesics do have the critical lower limit of $r_c=6M$ : this is a consequence of the presence of a black hole rather than just the Newtonian singularity.

The interactive radial motion is

While bound orbits exist in principle for any initial radial position larger than $r_c=6M$ , in practice the potential well becomes extremely flat for large values of this variable. Since this makes it difficult to provide a stable visualization in this region, and because the associated long periods are rather uninteresting, the upper extent of the radial variable is here capped at $r_0=100$ . At the lower end of its range the displayed value may differ slightly from that actually used in the evaluation in order to ensure the circular orbit at the minimum of the energy potential.

This behavior is now extended to three dimensions straighforwardly:

Given that the geodesics are generally not closed, while one could present multiple geodesics simultaneously, as was done in the presentations mentioned at the outset, this would get confusing rather quickly. Instead of doing that, note some variations in initial parameters that produce various sets of related geodesics.

Geodesics are rotated about the z-axis by adjusting ${\phi }_0$ , as previously noted. They are rotated about the apastron, the initial point furthest from the black hole, by decreasing the partial angular momentum $L_z$ from its maximum value for the given initial angle. They can be rotated about the periastron, the point of closest approach to the black hole, by varying ${\theta }_0$ while keeping the partial angular momentum at its maximum for that initial angle.

Having spatial visualizations of geodesics, one can easily show temporal variations along the vertical axis. Given that there is an extra spatial coordinate, a projection of some sort must be employed. The simplest is an orthographic projection in which the z-coordinate is ignored:

Since this interactive graphic retains all of the inputs from the previous graphic, one can in principle manipulate it to display temporal evolution in any desired invariant plane from any angle.

Possible additions for the future include animations of geodesic motion, which would require updates to the visualization library.

This presentation continues with the Kerr metric here.

Uploaded 2026.10.12 analyticphysics.com